﻿#pragma once
#pragma once
#include <stdio.h>
#include <xmmintrin.h>
#include <mmintrin.h>
#include <immintrin.h>
#include <smmintrin.h>
#include <intrin.h>  
namespace CoreCalc{
	class powerApproximation {

	public:
		powerApproximation() {}
		~powerApproximation() {}

	private:

	public:
		/* ————————————————————–  */
		/* (C)Copyright 2001,2008,                                         */
		/* International Business Machines Corporation,                    */
		/* Sony Computer Entertainment, Incorporated,                      */
		/* Toshiba Corporation,                                            */
		/*                                                                 */
		/* All Rights Reserved.                                            */
		/*                                                                 */
		/* Redistribution and use in source and binary forms, with or      */
		/* without modification, are permitted provided that the           */
		/* following conditions are met:                                   */
		/*                                                                 */
		/* - Redistributions of source code must retain the above copyright*/
		/*   notice, this list of conditions and the following disclaimer. */
		/*                                                                 */
		/* - Redistributions in binary form must reproduce the above       */
		/*   copyright notice, this list of conditions and the following   */
		/*   disclaimer in the documentation and/or other materials        */
		/*   provided with the distribution.                               */
		/*                                                                 */
		/* - Neither the name of IBM Corporation nor the names of its      */
		/*   contributors may be used to endorse or promote products       */
		/*   derived from this software without specific prior written     */
		/*   permission.                                                   */
		/*                                                                 */
		/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND          */
		/* CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES,     */
		/* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF        */
		/* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE        */
		/* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR            */
		/* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,    */
		/* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT    */
		/* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;    */
		/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)        */
		/* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN       */
		/* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR    */
		/* OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,  */
		/* EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              */
		/* ————————————————————–  */

#define madv(x, y, z)   _mm_add_ps(_mm_mul_ps(x, y), z)  
#define madd(x, y, c1)  _mm_add_ps(_mm_mul_ps(x, y), _mm_set_ps1(c1))  
#define madc(x, c0, c1) _mm_add_ps(_mm_mul_ps(x, _mm_set_ps1(c0)), _mm_set_ps1(c1))  

		/* 2^x: define if you want to raise 2 to positive values
		* x^y: define if you want to raise values above 1
		*
		* disables bias correction and overflow-check
		*/
#undef  EXP_POSITIVE_X  

		/* 2^x: define if you want to raise 2 to positive values
		* x^y: define if you want to raise values [0,1] by more than 258
		*
		* replaces underflow-check by max()
		*/
#undef  EXP_HUGE_X  

		/*
		* FUNCTION
		*    __m128 _exp2f4(__m128 x)
		*
		* DESCRIPTION
		*    The _exp2f4 function computes 2 raised to the input vector x.
		*      Computation is performed by observing the 2^(a+b) = 2^a * 2^b.
		*    We decompose x into a and b (above) by letting.
		*    a = ceil(x), b = x - a;
		*
		*    2^a is easilty computed by placing a into the exponent
		*    or a floating point number whose mantissa is all zeros.
		*
		*    2^b is computed using the following polynomial approximation.
		*    (C. Hastings, Jr, 1955).
		*
		*                __7__
		*          \
		*           \
		*    2^(-x) =   /     Ci*x^i
		*                /____
		*                 i=1
		*
		*    for x in the range 0.0 to 1.0
		*
		*    C0 =  1.0
		*    C1 = -0.9999999995
		*    C2 =  0.4999999206
		*    C3 = -0.1666653019
		*    C4 =  0.0416573475
		*    C5 = -0.0083013598
		*    C6 =  0.0013298820
		*    C7 = -0.0001413161
		*
		*/
		static __m128 exp2f4(__m128 &x) {
			__m128 SM_LN2 = _mm_set_ps1(0.69314718055994530942f);   /* log_e 2 */
			__m128 frac, frac2, frac4;
			__m128 result;
			__m128 hi, lo;

			__m128 exp_int, exp_frac;
			__m128i ix;

			/* for pow([0.0,1.0],[0.0,X])
			* -> log2([0.0,1.0]) * [0.0,X]
			* -> [-INF,-0.0] * [0.0,X]
			*
			* -INF  -> bias = 0×00000000, ix = -0×80000000, underflow = f, overflow = f
			* -5.00 -> bias = 0×00000000, ix = -5         , underflow = t, overflow = f
			* -0.20 -> bias = 0×00000000, ix =  0         , underflow = t, overflow = f
			* -0.00 -> bias = 0×00000000, ix =  0         , underflow = t, overflow = f
			*  0.00 -> bias = 0×3F7FFFFF, ix =  0         , underflow = t, overflow = f
			*
			* maximum “log2(smallest value not 0)” for 1/256 is -8
			* maximum “log2(smallest value not 0)” for 1/65536 is -16
			* maximum “log2(0)” for 0 is -127
			*/

			/* Break in the input x into two parts ceil(x), x - ceil(x).
			*/
#if   defined(EXP_POSITIVE_X)  
			__m128 bias;

			bias = _mm_cmpge_ps(x, _mm_set_ps1(0.0f));
			bias = _mm_and_ps(bias, _mm_castsi128_ps(_mm_set1_epi32(0×3F7FFFFF)));
			ix = _mm_cvtps_epi32(_mm_add_ps(x, bias));
#else  
			ix = _mm_cvtps_epi32(x);
#endif  

			frac = _mm_sub_ps(_mm_cvtepi32_ps(ix), x);
			frac = _mm_mul_ps(frac, SM_LN2);

#if   defined(EXP_POSITIVE_X)  
			__m128i underflow, overflow;
			__m128i exponent;

			underflow = _mm_cmpgt_epi32(ix, _mm_set1_epi32(-128));
			overflow = _mm_srli_epi32(_mm_cmpgt_epi32(ix, _mm_set1_epi32(128)), 1);
			exponent = _mm_add_epi32(ix, _mm_set1_epi32(127);
			exponent = _mm_slli_epi32(exponent, 23);
			exponent = _mm_and_si128(exponent, underflow);
#elif   defined(EXP_HUGE_X)  
			__m128i underflow;
			__m128i exponent;

			underflow = _mm_cmpgt_epi32(ix, _mm_set1_epi32(-128));
			exponent = _mm_add_epi32(ix, _mm_set1_epi32(127);
			exponent = _mm_slli_epi32(exponent, 23);
			exponent = _mm_and_si128(exponent, underflow);
#else  
			__m128i exponent;

			/* log(0) -> -127 * 258 -> -32766 -> -127
			* log(0) -> -127 * 259 ->  32643 -> 32643
			*
			* maximum is 0^258
			*/
			exponent = _mm_max_epi16(ix, _mm_set1_epi32(-127));
			exponent = _mm_slli_epi32(exponent, 23);
			exponent = _mm_add_epi32(exponent, _mm_set1_epi32(127 << 23));
#endif  

			exp_int = _mm_castsi128_ps(exponent);

			/* Instruction counts can be reduced if the polynomial was
			* computed entirely from nested (dependent) fma's. However,
			* to reduce the number of pipeline stalls, the polygon is evaluated
			* in two halves (hi amd lo).
			*/
			hi = madc(frac, -0.0001413161f, 0.0013298820f);
			lo = madc(frac, -0.1666653019f, 0.4999999206f);

			frac2 = _mm_mul_ps(frac, frac);

			hi = madd(frac, hi, -0.0083013598f);
			lo = madd(frac, lo, -0.9999999995f);

			frac4 = _mm_mul_ps(frac2, frac2);

			hi = madd(frac, hi, 0.0416573475f);
			lo = madd(frac, lo, 1.0f);

			exp_frac = madv(frac4, hi, lo);

			result = _mm_mul_ps(exp_frac, exp_int);

#ifdef  EXP_POSITIVE_X  
			/* Handle overflow (+INF) */
			result = _mm_or_ps(result, _mm_castsi128_ps(overflow));
#endif  

			return (result);
		}

		/*
		* FUNCTION
		*    vector float _log2f4(vector float x)
		*
		* DESCRIPTION
		*    The _log2f4 function computes log (base 2) on a vector if inputs
		*      values x. The _log2f4 function is approximated as a polynomial of
		*      order 8 (C. Hastings, Jr, 1955).
		*
		*                   __8__
		*             \
		*              \
		*    log2f(1+x) =  /     Ci*x^i
		*                   /____
		*                    i=1
		*
		*    for x in the range 0.0 to 1.0
		*
		*    C1 =  1.4426898816672
		*    C2 = -0.72116591947498
		*    C3 =  0.47868480909345
		*    C4 = -0.34730547155299
		*    C5 =  0.24187369696082
		*    C6 = -0.13753123777116
		*    C7 =  0.052064690894143
		*    C8 = -0.0093104962134977
		*
		*    This function assumes that x is a non-zero positive value.
		*
		*/
		static __m128 log2f4(__m128 &x) {
			__m128i exponent, uxponent;
			__m128 result;
			__m128 x2, x4;
			__m128 hi, lo;
			__m128 y;

			/* for pow([0.0,1.0],[0.0,X])
			* -> log2([0.0,1.0])
			* -> [-INF,-0.0]
			*
			* 1.00 -> exponent =    0, y = 0.0
			* 0.25 -> exponent =   -2, y = 0.0
			* 0.01 -> exponent =   -7, y = 0.28
			* 0.00 -> exponent = -127, y = 0.0
			*/

			/* Extract the exponent from the input X.
			*/
			exponent = _mm_and_si128(_mm_castps_si128(x), _mm_set1_epi32(0×7f800000));
			uxponent = _mm_sub_epi32(exponent, _mm_set1_epi32(127 << 23));
			exponent = _mm_srai_epi32(uxponent, 23);

			/* Compute the remainder after removing the exponent.
			*/
			y = _mm_castsi128_ps(_mm_sub_epi32(_mm_castps_si128(x), uxponent));

			/* Calculate the log2 of the remainder using the polynomial
			* approximation.
			*/
			y = _mm_sub_ps(y, _mm_set_ps1(1.0f));

			/* Instruction counts can be reduced if the polynomial was
			* computed entirely from nested (dependent) fma's. However,
			* to reduce the number of pipeline stalls, the polygon is evaluated
			* in two halves (hi amd lo).
			*/
			hi = madc(y, -0.0093104962134977f, 0.052064690894143f);
			lo = madc(y, 0.47868480909345f, -0.72116591947498f);

			x2 = _mm_mul_ps(y, y);

			hi = madd(y, hi, -0.13753123777116f);
			lo = madd(y, lo, 1.4426898816672f);

			x4 = _mm_mul_ps(x2, x2);

			hi = madd(y, hi, 0.24187369696082f);
			lo = _mm_mul_ps(y, lo);

			hi = madd(y, hi, -0.34730547155299f);

			result = madv(x4, hi, lo);

			/* Add the exponent back into the result.
			*/
			result = _mm_add_ps(result, _mm_cvtepi32_ps(exponent));

			return (result);
		}

		static __m128 powf4(__m128 &x, __m128 &y) {
			/* for pow([0.0,1.0],[0.0,X])
			* -> log2([0.0,1.0]) * [0.0,X]
			* -> log2(0,00390625) * 2.2    -> -8
			* -> log2(0,00390625) * 0.45   -> -8
			* -> [-INF,-0.0] * [0.0,X]
			* -> [-127,-0.0] * [0.0,X]
			*/
			return exp2f4(_mm_mul_ps(log2f4(x), y));
		}

		static __m128 powf4(__m128 &x, float n) {
			__m128 y; y = _mm_set_ps1(n);

			return powf4(x, y);
		}

		static __m128 powf4(__m128 &x, long n) {
			__m128 y; y = _mm_cvt_si2ss(y, n);

			return powf4(x, y.m128_f32[0]);
		}

#undef  madv  
#undef  madd  
#undef  madc  

	};


	ref class Funcs
	{
	public:

		Funcs()
		{
		}
	};

}